

namespace edu.stanford.hci.ideas
{

    /// <summary>
    /// Summary description for TransPanel.
    /// </summary>
    public class FFT
    {
        double[] outputData;
        public void SWAP (ref double x, ref double y) {
            double temp;
            temp = y;
            y = x;
            x = temp;
        }

        public double GetAmplitude(int index)
        {
            return System.Math.Sqrt(outputData[2*index] * outputData[2*index] + outputData[2*index + 1] * outputData[2*index + 1]);
        }

        public ulong ComplexFFT(double[] data, ulong number_of_samples, ulong sample_rate, int sign)
{

    outputData = new double[2 * sample_rate];
    for (int counter = 0; counter < 2 * (int) sample_rate; counter++)
    {
        outputData[counter] = 0;
    }
    ulong fundamental_frequency = 0;
        //variables for the fft
        ulong n,mmax,m,j,istep,i;
        double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
                
        //the complex array is real+complex so the array
    //as a size n = 2* number of complex samples
    //real part is the data[index] and
    //the complex part is the data[index+1]
 
        //new complex array of size n=2*sample_rate
      
      
        
        //put the real array in a complex array 
     //put the real array in a complex array
        //the complex part is filled with 0's
        //the remaining outputData with no data is filled with 0's
        for(n=0; n<sample_rate;n++)
        {
                if(n<number_of_samples)
                        outputData[2*n]=data[n];
                else
                        outputData[2*n]=0;
                outputData[2*n+1]=0;
        }
        
        //binary inversion (note that the indexes
    //start from 0 witch means that the
    //real part of the complex is on the even-indexes
    //and the complex part is on the odd-indexes)
        n=sample_rate << 1;
        j=0;
       for (i=0;i<n/2;i+=2) {
                if (j > i) {
                        SWAP(ref outputData[j],ref outputData[i]);
                        SWAP(ref outputData[j+1],ref outputData[i+1]);
                        if((j/2)<(n/4)){
                                SWAP(ref outputData[(n-(i+2))],ref outputData[(n-(j+2))]);
                                SWAP(ref outputData[(n-(i+2))+1],ref outputData[(n-(j+2))+1]);
                        }
                }
                m=n >> 1;
                while (m >= 2 && j >= m) {
                        j -= m;
                        m >>= 1;
                }
                j += m;
        }
        //end of the bit-reversed order algorithm
   //Danielson-Lanzcos routine
        mmax=2;
        while (n > mmax) {
                istep=mmax << 1;
                theta=sign*(2*System.Math.PI/mmax);
                wtemp=System.Math.Sin(0.5*theta);
                wpr = -2.0*wtemp*wtemp;
                wpi=System.Math.Sin(theta);
                wr=1.0;  
                wi=0.0;
                for (m=1;m<mmax;m+=2) {
                        for (i=m;i<=n;i+=istep) {
                                j=i+mmax;
                                tempr=wr*outputData[j-1]-wi*outputData[j];
                                tempi=wr*outputData[j]+wi*outputData[j-1];
                                outputData[j-1]=outputData[i-1]-tempr;
                         outputData[j]=outputData[i]-tempi;
                                outputData[i-1] += tempr;
                                outputData[i] += tempi;
                        } 
                        wr=(wtemp=wr)*wpr-wi*wpi+wr;
                        wi=wi*wpr+wtemp*wpi+wi;
                }
                mmax=istep;
        }
        //end of the algorithm
                
        //determine the fundamental frequency
        //look for the maximum absolute value in the complex array
        fundamental_frequency=0;
        for(i=2; i<=sample_rate; i+=2)
        {
       if((System.Math.Pow(outputData[i],2)+System.Math.Pow(outputData[i+1],2))>(System.Math.Pow(outputData[fundamental_frequency],2)+System.Math.Pow(outputData[fundamental_frequency+1],2))){

                        fundamental_frequency=i;
                }
        }
                
        //since the array of complex has the format [real][complex]=>[absolute $
        //the maximum absolute value must be ajusted to half
        fundamental_frequency=(ulong)System.Math.Floor((float)fundamental_frequency/2);
        return fundamental_frequency;
}
    }
}
